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Abstract 

We introduce a new method to perform preliminary orbit determination 
for space debris on low Earth orbits (LEO). This method works with tracks 
of radar observations: each track is composed by n > 4 topocentric position 
vectors per pass of the satellite, taken at very short time intervals. We as¬ 
sume very accurate values for the range p, while the angular positions (i.e. 
the line of sight, given by the pointing of the antenna) are less accurate. We 
wish to correct the errors in the angular positions already in the computa¬ 
tion of a preliminary orbit. With the information contained in a pair of radar 
tracks, using the laws of the two-body dynamics, we can write 8 equations 
in 8 unknowns. The unknowns arc the components of the topocentric ve¬ 
locity orthogonal to the line of sight at the two mean epochs of the tracks, 
and the corrections A to be applied to the angular positions. We take advan¬ 
tage of the fact that the components of A are typically small. We show the 
results of some tests, performed with simulated observations, and compare 
this algorithm with Gibbs’ method and the Keplerian integrals method. 


1 Introduction 

We investigate the preliminary orbit determination problem for a satellite of the 
Earth using radar observations collected by an instrument with given technical 
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specifications, and with a fixed observation scheduling. Assume we collect the 
following data for the observed object: 

( tj,pj,aj,6j ), j - 1...4 (1) 

where the triples (py, a,, dy) represent topocentric spherical coordinates of the ob¬ 
ject at epochs tj. Typically cry, dy are the values of right ascension and declination. 
We shall call radar track the set of observations in (1). 

The following assumptions will be made on the data composing the tracks. 
The time difference tj+\ - tj between consecutive observations is At = 10 s. The 
range data pj are very precise: the statistical error in the range is given by its RMS 
o~p, which is 10 m. On the other hand we assume that the angles ay, dy are not 
precisely determined: their RMS cr a , cr s are supposed to be 0.2 degrees. 

Given a radar track we can compute by interpolation the following data: 

(t,d,8,p,p,p). (2) 

Here t , a and 8 are the mean values of the epoch and the angles, and p, p, p are 
the values of a function pit) and its derivatives at t -t, where pit) is given by a 
quadratic fit with the (ty,py) data. 

For low Earth orbits (LEO) these assumptions imply that the interpolated val¬ 
ues of dr, 8 are very badly accurate, to the point that their value can be of the same 
order of the errors, therefore they are practically undetermined. 

By the above considerations, given a vector (2) obtained by a radar track, and 
using spherical coordinates and velocities 

(p, a,8,p, dr, 8) 

to describe the orbit, we can consider as unknowns the quantities (Aa, Ad, dr, 8), 
with 

a = a + Aa, 8 = 8 + Ad, 

where A a, A 8 are small deviations from the mean values d, 8. 

To search for the values of the unknowns we need to use additional data: we 
can try to use the data of 2 radar tracks, together with a dynamical model, to 
compute one or more preliminary orbits. This is a linkage problem, see [5]. 

In this paper we propose a new algorithm for the linkage, which takes ad¬ 
vantage of the smallness of A a, Ad, that we call infinitesimal angles. We write 
the equations for preliminary orbits by using the 5 algebraic integrals of Kepler’s 


2 



problem, Lambert’s equation for elliptic motion (see Section 12) and the projec¬ 
tion of the equations of motion along the line of sight. 

Moreover, we perform some tests to compare this method with Gibbs’ method, 
using only one radar track, and with the Keplerian integrals method, which solves 
a linkage problem using (a, 6,p,p) at two mean epochs (see [9], [2], [3]). 

The paper is organized as follows. First we introduce some notation and recall 
the basic results on Kepler’s motion which are relevant for this work (see Sec¬ 
tions 2, 3, 4). The equations for the linkage problem, see (7), are presented in 
Section 5, and in Sections 6, 7 we show two different ways to compute solutions 
of (7). In Section 9 we present the results of some numerical tests, including a 
comparison with the already known methods recalled in Section 8. Finally, in 
Section 12, we recall the proof of Lambert’s theorem for elliptic orbits and give 
a geometrical interpretation of the results. Moreover, we show a method to cor¬ 
rect the observations of a radar track so that they correspond to points in the same 
plane. 


2 The equations of motion 

Let us denote by e ; ’ the unit vector corresponding to the line of sight, and by q 
the geocentric position of the observer. Then the position of the observed body is 
r = q + pe p , where p is the range. Using the right ascension a and the declination 
5 as coordinates we have 



3 



where r/ = Vd 2 cos 2 6 + S 2 is the proper motion and k = • n. For later use we 

introduce the notation 

7C = (r + -^-r) • = p - prj 1 + q • e^’ + -^-(r • e p ). 

V | r |J / | r |J 


3 The two-body integrals 


We write below (see also [3]) the expressions of the first integrals of Kepler’s 
problem, i.e. the angular momentum c, the energy & and the Laplace-Lenz vector 
L, in the variables p, a, 5,p, f f with 

^ = pa cos 6, £ = pi). (4) 


We have 


where 


with 


and 


c = + B4 + C, 



pL(p,p) = r x c -{*y- = (|r| 2 - n)r - (r • r)r, 
|r| v |r| / 

A = rxe", B = r x e' s , C = rxq + J 6qxe p , 

a 1 ^ , d& 

e “ cos 6 da’ 6 “ 06' 


r = <fe" + Cf + (peP + q), 

|r| 2 = f 2 + if + 2q • + 2q • e 5 £ + | pe p + q| 2 , 

r r = q e a % + q • e 6 £ + (peP + q) • r. 


We introduce the notation 

q a = q • e“, q s = q • e 6 , if = q • e“, if - q • e 5 . 
Note that f 2 + if - p 2 r/ 2 . 
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4 Lambert’s equation 


Lambert’s theorem for elliptic motion gives the following relation for the orbital 
elements of a body on a Keplerian orbit at epochs t\,t 2 . 

n(t 2 - h) = ft - y - (sinft - siny) + 2 kn. (5) 


Here k 6 N is the number of revolutions in the time interval [t\, t 2 ], n = n(a) is the 
mean motion, where a = -/r/(2fi) (the energy is the same at the two epochs), and 
the angles ft. y are defined by 


, ft r, + r 2 + d 


snr - = 


4 a 



r \ + r 2 — d 
4a 


( 6 ) 


and 

0 < ft - y < 2n, 


with /-|, r 2 the distances from the center of force, and d the length of the chord 
joining the two positions of the body at epochs t\,t 2 . For a fixed number of revo¬ 
lutions we have 4 different choices for the pairs (J3,y ), see Section 12 and [1] for 
the details. 


5 Linkage 

We wish to link two sets of radar data (2), with mean epochs tj, i = 1,2, and 
compute one or more preliminary orbits. In the following we use labels 1, 2 for 
the quantities introduced in Sections 2, 3, 4 according to the epoch. 

Let us denote by X the expression defining Lambert’s equation. More pre¬ 
cisely, X = 0 is one of the possible cases occurring in (5), see Section 12. More¬ 
over, let us define \ 2 = x q 2 . We consider the system 

(Ci — c 2 ,fii — S 2 , f K\, c K 2 ,(L\ — L 2 ) • v 2 ,X) = 0 (7) 

of 8 equations in the 8 unknowns (X, A), with 

X = (£i, £i, O 2 X 2 X A = (Acri, Adi, Aq; 2 , Ad 2 ). 

Note that the unknowns are divided into 2 sets so that A is the vector of infinites¬ 
imal angles. To solve system (7) we first compute X as function of A using 4 of 
these equations, then we substitute X(A) into the remaining equations and search 
for solutions of the resulting nonlinear system by applying Newton-Raphson’s 
method. Taking advantage of the assumed smallness of the solutions A, we can 
use A = 0 as starting guess. 
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6 Computing X(A) 

We describe below two methods to compute X as function of A using some of 
the equations of system (7). One approach uses linear equations, see Section 6.1, 
while the equations for the other are quadratic, see Section 6.2. 

6.1 Linear equations 

Substituting 2<5| + p{K\ - 2£? - PiKi in place of &\ - S 2 in (7) we obtain an 
equivalent system and the equation 


2fi] + p\K\ — 2£ 2 + Pi^2 


( 8 ) 


is linear in the variables X = (£ 1; f 2 , fi)- 

Using equation (8) and the conservation of the angular momentum we obtain a 
linear system in the variables X: 


Here 


M = 


mx = v. 


An 

B 11 

-A 21 

~B 2 \ 

Ai 2 

7?12 

-a 22 

~b 22 

A 13 

5i3 

-A?3 

-B 23 




-q 2 


(9) 


where A, P ] are the components of A, , B/, and q a { - q ( --e“, if - q for z = 1,2. 
Moreover 


V - (C 21 - Cn,C 22 - Ci 2 ,C 2 3 - C l3 ,D 2 - Df ) T , 


where Cq are the components of C ; and 


Di = 


UpWj + I P«< + q,| 2 ) 


P_ 

|r,l 


with rf expressed as function of (Aar,-, Adf) by using the equations TC, = 0, i = 1,2, 
that is using relation 

if = ~(p + q • (r • o) 

at the 2 epochs t \, t 2 . 
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We can write X as function of A by solving system (9). 

Let us call Mhj the components of M, and V/, the components of V. The solutions 
of (9) are given by 


\M\ ’ 


where M/ ( has components 


\Mii\ 
\M\ ’ 


f = 1,2 



yVl/j j if k ± j 

14 if k = j 


and \M\, \Mk\ represent the determinants of M, Al/,. 


( 10 ) 


6.2 Quadratic equations 

The orbits at epochs t \, 4, computed with the solution X of system (9), do not 
necessarily share the same energy <5. This can produce some problems in the 
linear algorithm described above, especially when solving Lambert’s equation, 
where the right-hand sides of ( 6 ) may become greater that 1 during the iterations 
of Newton-Raphson’s method. We can force the orbits to share the same energy 
by solving the first 4 equations in (7), that are quadratic equations in the variable 
X. 

By introducing the vector 

Y = (6,6,&), 

we can write the conservation of the angular momentum as the linear system 


Here 


and 


where 



NY 

= w. 


(ID 


' A n 

flu 

-A 21 


N = 

An 

Bn 

-A 22 

(12) 


_ A ,3 

Bn 

-A 23 


W 

= £W (1) + W (0) , 



- (#21, ^22, B 2 3 ) T , 

= (C21 - Cn, C22 - C\ 2 , C23 - C\i) T . 


w (1) 

w (0) 
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We solve system (11). Let us call N h j the components of N and W h ,Wf ] ,W^ ] the 
components of W,W (0) ,W (1) . The solutions of (11) are functions of 6, A, and are 
given by 

~ imi ~ m ~ m 

^ |W| ’ 41 \N\ ’ 4:2 |JV| ’ 

where A4 has components 

sj(k) = / N h j if k±j 
hj \W h if k = j ■ 

From the conservation of energy we can find 6 as function of A. We write 

F 2 £ + F l £ 2 + F 0 = 0, (13) 

with 

f 2 = j^(i<¥ + i<i 2 - i<i 2 ) -1 

F \ = ^(l<||<| + IW^IITVf | - IJVfllA/f |) + 

+ + <A\n?\ - q a i l< } l - q\m 

F ° = jTTpd^ri 2 + K 0) l 2 - l^fi 2 ) 

+ X)i — X>2, 

where N { k c \ k = 1,2,3, £ = 0,1, has components 


and 


Therefore we have 


= 

hj 


Nhj if k±j 
wf if k = j ’ 


$«■ = 2 A - ppit 


i = 1,2. 


6(A) = 6 (6(A), A), 
6(A) = 6 (6(A), A), 
6(A) = 6(6(A), A), 


where 6(A) is a solution of (13). Note that we can have up to two acceptable 
expressions for X(A). 
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7 Computing A 

We introduce the vector 


G — CK\,% 2 , (Li — L 2 ) • v 2 , X). 


Remark 1 To select the relevant expressions of X H’<? need to guess the value of 
k in (5). We can do this by assuming A = 0 and computing the possible orbits 
according to the linear or quadratic equations for X(A). In both cases we obtain 
two possible values for the number of revolutions k: with the linear equations we 
can obtain two different values of k at the two epochs t\, h; with the quadratic 
equations we may obtain two orbits with different k at the same epoch, say f, but 
from conservation of energy we obtain the same values at t 2 . 

By substituting the possible expressions of X(A), coming from either the linear or 
the quadratic equations, we obtain the reduced system 

£(A) = G(X(A),A) = 0. (14) 


Since the unknowns in A are small, we can try to apply Newton-Raphson’s method 
with A = 0 as starting guess. Thus we try to compute an approximation for A by 
the iterative formula 

A, +1 = A, - [^|(A A )f £(A a ), A 0 = 0. (15) 

Equations (15) are linear, and are defined by (14) and by the Jacobian matrix 


d -f (A*) = f(X„A„)f(A*) + ^(X*, A*), 
aA oX o A o A 


with X k = X(A a ). 

Remark 2 Note that at each iteration the number of solutions can be doubled, 
but if we impose the value of A a+1 to be close to \ k then we can usually avoid 
bifurcations. 

The computation of the Jacobian matrix || is described below, enhancing the 
differences between the linear and the quadratic case. 
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7.1 The derivatives 


d<K x 

~dX 

d'K, 

~dX 


-( 6 ,£, 0 , 0 ) 

Pi 

--(0,o,&, 6) 

P2 


We observe that 


Thus we have 


L 2 • v 2 = --(r 2 • r 2 )(r 2 • v 2 ). 
P 


——[(Li - L 2 ) • v 2 ] 


— [(Li - L 2 ) • v 2 ] 


— [(Li - L?) • v 2 ] 

[(Li - L 2 ) • v 2 ] 

dCi 


■■ -(^1 + qi • e")(r! • v 2 ) 

P 

--[(qi • e“)(r! • v 2 ) + (r, • rOCe? • v 2 )] 

n 

■■ -(£ + qi • ef)(r! • v 2 ) 

P 

--[(qi • ej)(ri • v 2 ) + (i-i • ri)(ej • v 2 )] 

JJ, 

■ - [(q2 • e 2 )(r 2 • v 2 ) + (r 2 • r 2 )(e 2 • v 2 )] 

/i 

: -[(q 2 • e^)(r 2 • v 2 ) + (r 2 • r 2 )(e^ • v 2 )] 
P 


For Lambert’s equation, the derivatives are given by 


d£ dn _ _ d(fi - sin /3) d(y - sin y) 

dX = dX {t '~ t2 ’ + dX dX ’ 

dn 3 r— — 5(2fii) 

PX - 


d(/3-sin f3) 

dX 

d(y - sin y) 

dX 


„ n ,d[3 „ / r + dY + 


„ Jy „ r - ar_ 
(1_cosr) 5x = 2 Vr^T:5x’ 
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with 


r+ 

r 




r\ + r 2 + d 
2[x 


r i + r 2 - d 

2[i 


In the expression for we use the energy t)\ at epoch t\. We could as well choose 
£ 2 at epoch t 2 . this choice is arbitrary in the linear case, in fact computing X(A) 
with the linear algorithm, we generally have £i(£i(A),£i(A)) ± £2 (£2 (A), £2 (A)). 
Since r 1? r 2 , d do not depend on X, we have 

3T + r { + r 2 + d d&\ 

~dX = 2 ~n dX’ 

dr r\ + r 2 — d d£>\ 

dX = 2/u ax’ 


with 


dS 1 
dX 


- ((£1 + Qi • e“), (("1 + ffi ' e i)’ 0,0). 


7.2 The derivatives H 


To compute derivatives with respect to A we use as intermediate variables the unit 

vectors ef 3 -, e". e' 5 , / = 1,2. To this aim we introduce the vector 

j j j J 





f e 1 ' \ 


' Ei ~ 

e 2 


j 

E = 

, where E j = 

pQ! 

y 

e 5 


Its derivatives with respect to A are given by 


dE 

dA 


8E l 

d(ai,Si ) 

0 


0 

8E 2 

8(0:2,62) - 


where 


gg j_ 

dips j, Sj ) 


cos 8 jt a . e 6 ; 

J J J 

e| 0 

-sind.e" -e p ; 

J J J 


and e| = -(cos aj, sin aj, 0) T . 

Moreover, we need to compute We describe the different procedures for 
the linear and quadratic methods. 
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7.2.1 The derivatives ||, linear case 

Using (10), we only need to compute 

56 = 1 d|AU/,-il |M 2 /,-il d\M\ 

dE \M\ dE |M| 2 dE ’ 

56 = 1 d\M 2h \ \M lh \d\M\ 

dE \M\ dE |M| 2 dE ' 

We take advantage of the following relation, valid for any matrix A of order n with 
coefficients a/j depending on a variable x: 


d_ 

dx 


|A| = I B h 

h= 1 


where B h has coefficients bf), with 



dx '' 


h F j 
h = j 


7.2.2 The derivatives ||, quadratic case 


From the implicit function theorem applied to equation (13) we obtain 


56 _ _ r 

dE ~ I- 


1 


(IF2^2 + F 1 )\ dE 


,dF^ 2 dF i dF oy | 

(*N + *r& + ~) 


ft=4°(E) 


Let us define 


6 (E) = 6(6( E),E), 
6(E) = 6(6(E),E), 
6(E) = 6(6(E),E). 


We have 


56 

- 

r 5|MI 

56 

d|Mk 

IMP 

5|M| 

dE 

IMP 

1 56 

dE 

dE ' 

IMP 

dE 


- 

r5|M 2 | 

56 

5|M 2 K 

m 

5|M| 

dE 

IMP 

v 56 

dE 

+ dE ' 

IMP 

dE 

56 

- 

<5|Ad 3 | 

56 

d|7V 3 K 

|M 3 | 

5|M| 

dE 

|M| 

1 56 

dE 

+ dE ' 

IMP 

dE 
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7.3 The derivatives ^ 

As in Section 7.2, we compute the derivatives of Q with respect to E and multiply 
the result by We have 


and 


PTC 




j _ .. 

- q/ + Pi„ 13 


l j i 


1 - 3 pj 


( r rO N 


7=1,2 


PTC/ _ PTC/ (77C, _ PTC, 

PeJ de 5 j ’ PE 2 PEi 


^p[(Li - L 2 ) • v 2 ] = -^[(r! • v 2 )| 2 p!qi+/rpi^J + 

+ (|ri | 2 - j^j)pi v 2 - (*1 • v 2 )(piqi +piqi) - (rj • ri)piv 2 ], 
P £\ 

— [(Lj -L 2 )-v 2 ] = — [2(rj-v 2 )qi-(r! • v 2 )qi-(rj-rOvs], 

d / 

— [(Lj -L 2 ) • v 2 ] = —[2(r! • v 2 )q!-(rj • v 2 )q!-(rj-rOv^, 

at " p 


P 

K 


[(Li - L 2 ) • v 2 ] - -Lj x q 2 + 


+ -[(P2q2 +P2q2)(f2 • V 2 ) + (r 2 • r 2 )q 2 x q 2 l, 
pL J 


A[( L ,-L 2 ).v 2 ] 


-[&(r 2 • v 2 ) + ^ 2 (r 2 • r 2 )]q 2 , 
P 

5 1 

— [(Li - L 2 ) • v 2 ] = -[£(r 2 • v 2 ) - £ 2 (r 2 • r 2 )]q 2 . 


<9e*2 P 

For Lambert’s equation we have 


P£ dn _ _ <9(jS - sin /?) P(y - sin y) 

PE = PE Ul " h ’ + PE PE ' 


with 


rWg~ d(2<Si) ^ 

PE! 2p V 1 PE! ’ PE 2 


= 0, 
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d(/3 - sin/3) o 

/ r + 0T + 

/IE 2 A 

1 1 - r + oe 

d(y - sin y) o 

"I 

dE “ 2 A 

/1 - r_ oe 


Moreover 

dr + _ 2 Gi.dn dd Y + 0(260 

dE[ ~ ~^lu { dE l + dE l )+ 26 { OE } ’ 

5r + 26 1 dr 2 dd 

SE 2 = ~^p { dE 2 + dE 2 

dr_ _ 2 Gi.dn dd T_ 0(260 

dE { ~ ~^u i 0 E l ~dE l )+ 26 l OEi ’ 

5r_ 261 Or 2 dd 

dE 2 = ~^p { dE 2 ~dE 2 


with 


9(260 . . _ qi . or ■ \ 

= (2p 1 q l + 2ppi—, 2^ r ] q 1 , 2frq0, 


<9E, 

£ = (—■«■«>■ 
oEj r, 

- r 2’°’°)’ 

£/Ei u 


^ = (^,0,0), 

oE 2 r 2 

= ^y(q 2 -r 1 , 0 , 0 ). 

<yE? a 


8 Alternative known methods 


We recall below two already known methods that can be used in place of the 
algorithm described in Sections 5, 6, 7, with the available data. An important 
difference is that these methods do not provide corrections to the angles a, 8. 


8.1 Gibbs’ method 

From three position vectors of an observed body at the same pass we can compute 
an orbit using Gibbs’ method, see [4, Chap. 8]. We recall below the formulas of 
this method. Given the position vectors r j,j = 1,2,3 at times tj, Gibbs’ method 
gives 

r 2 = -d\ ri + d 2 r 2 + d 3 r 3 . 
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where 


dj 

G] 


Gj + Hj,y 


'32 


, G 3 


j= 1,2,3, 

f 

'21 


, G 2 = G i - G 3 , 


h\hih\ h\hih\ 

p t 32 /U , h 3 = ntn! 12 ,h 2 = h { - h 3 . 


Here t,j = 0 - tj, rj = |r ; -|. 


8.2 Keplerian integrals 

From two radar tracks, we can obtain by interpolation the values of ( a,6,p,p ) 
at epochs t \, t 2 . If we wish to determine the values of the unknowns a, 6, or 
equivalently of £ defined by (4), we can use the Keplerian integrals method, see 
[9], [2], [3]. This method uses the equations 

Cl = c 2 , &\ = & 2 , 

which can be explicitly solved, giving at most two solutions. 


9 Numerical tests 

We have performed some numerical tests with simulated objects, without the J 2 
effect, but adding errors to the observations. Here we describe the results of our 
tests only for one simulated object. We add to the data of the tracks a Gaussian 
error, with zero mean and standard deviations listed in Table 1. In particular we 
consider the case where we add no error to the range p (RMS (1) in the table). The 


RMS 

a 

5 

P 

(1) 

0.2 

0.2 

- 

(2) 

0.1 

0.1 

5 x 10^ 3 

(3) 

0.2 

0.2 

10" 2 


Table 1: RMS of the errors added to the radar tracks. 

data that we obtain by interpolation of two radar tracks are displayed in Table 2 
for the simulated object. The case labeled with RMS (1) is peculiar, in fact we 
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Epoch 

Data 

RMS (1) 

RMS (2) 

RMS (3) 


a (deg) 

51.17 

51.20 

51.18 


8 (deg) 

-5.47 

-5.44 

-5.47 

54127.155035 

p (km) 

1984.4 

1984.4 

1984.4 

(MJD) 

p (km/d) 

-73313 

-73268 

-73223 


p (km/d 2 ) 

116444362 

116534247 

116676527 


d (deg) 

264.30 

264.28 

264.30 


8 (deg) 

-66.77 

-66.79 

-66.77 

54127.582118 

p (km) 

1893.5 

1893.5 

1893.5 

(MJD) 

p (km/d) 

-323582 

-323712 

-323842 


p (km/d 2 ) 

123666885 

123168829 

122613669 


Table 2: Data interpolated from the radar tracks of the test object. 


interpolate the available values of a, 6 and we use the exact values of p,p,p, that 
we can compute from the given orbit. 

In Table 3 we show the orbits computed by the methods of Gibbs (G), by the Kep- 
lerian integrals (KI) and by the infinitesimal angles (InfAng), using the quadratic 
equations introduced in Section 6.2. For KI and InfAng we use the three data sets 
displayed in Table 2. Note that InfAng with RMS (1) is able to correct the errors 
in a, 8 and to recover the orbital elements of the known orbit. For RMS (2) and 
RMS (3), InfAng obtains a better value of the semimajor axis a, and slightly worse 
values of the other orbital elements, if compared with KI. To be consistent, for KI 
with RMS (1) we use the exact values of p,p. The results with Gibbs’ method are 
not very good. On the other hand it uses only part of the information contained in 
the data: here we use the three vectors (tj,pj, ay , 8j) of the first track at epochs tj, 
with j = 1,2,4. 

The infinitesimal angles method with linear equations shows some limitations, in 
a way that we could not compute an orbit for reliable values of the observational 
errors. 


10 Conclusions 

We have introduced a new method to compute preliminary orbits of space de¬ 
bris using radar observations. The comparison with already existing methods was 
performed for some test cases. Large scale tests should be done, to check the per¬ 
formance of the algorithm, possibly with real data. We plan to investigate the case 
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Known 

orbit 

G 

OD methods 

KI 

InfAng 

RMS 



8267.75 

7816.61 

7818.10 

(1) 

a 

7818.10 

8022.40 

7815.71 

7818.02 

(2) 



8267.91 

7814.80 

7818.09 

(3) 



0.005 

0.066 

0.066 

(1) 

e 

0.066 

0.037 

0.066 

0.065 

(2) 



0.005 

0.066 

0.065 

(3) 



70.17 

65.85 

65.81 

(1) 

I 

65.81 

68.02 

65.84 

65.61 

(2) 



70.17 

65.85 

65.34 

(3) 



217.55 

216.25 

216.25 

(1) 

D. 

216.25 

216.92 

216.25 

216.30 

(2) 



217.55 

216.25 

216.37 

(3) 



336.76 

357.00 

357.16 

(1) 

CO 

357.16 

356.16 

357.17 

357.43 

(2) 



335.93 

357.19 

357.53 

(3) 



219.54 

202.29 

202.08 

(1) 

t 

202.08 

201.46 

202.09 

201.72 

(2) 



220.38 

202.09 

201.53 

(3) 


Table 3: Orbital elements at epoch t\ = 54127.15 MJD computed by the tracks 
producing the three data sets of Table 2. Distances are expressed in km, angles in 
degrees. 

which includes the J 2 effect in the equations: this is essential to link radar tracks 
of LEO orbits after several revolutions. 
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12 Appendix 

12.1 Lambert’s equations for elliptic motion 

In this section, following [10] and [7], we summarize the steps to derive Lambert’s 
equation for elliptic motion under a Newtonian force and we give a geometric 
interpretation of the result. Indeed we obtain four distinct equations per number 
of revolutions of the observed body. Note that, dealing with radar observations 
of space debris, the time between two distinct arcs of observations usually covers 
several revolutions. 


Theorem 1 (Lambert, 1761) In the elliptic motion under the Newtonian gravita¬ 
tional attraction, the time At = h - t\ spent to describe any arc (without multiple 
revolutions) from the initial position P\ to the fined position Pi depends only on 
the semimajor axis a, on the sum r = r\ + r 2 of the two distances f| = \P\ - F\, 
r 2 = |P 2 ~ F\from the center of force F, and the length d of the chord joining P\ 
and Pi- More precisely we have 

nAt = fi - y - (sin fi - sin y), 

where n = n(a) is the mean motion, and the angles j3, y are defined by 


. 2 fi r + d 
Sm 2 = ~4a~’ 


. 2 y _r-d 
Sm 2 4a ’ 


0 < yS - y < 2 n. (16) 

Proof We can assume, without loss of generality, that the positions of the points 
Pi, Pi are defined by two values E\,E 2 of the eccentric anomalies such that 0 < 
Ei- E\ < 2n. 

The difference of Kepler’s equations at the two epochs gives 

nAt = E 2 - Pi - ^sinP^ - sin Pi), 


where e is the orbital eccentricity. From elementary geometrical relations we 
obtain 

r n,(-, E\ + Ei Ei-E { \ 

a \ 2 2 1 


d 0 . Ei- 

- = 2 sm-— 

a 2 


1 - e 2 cos 


2 Ei + Ei 
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It follows that 


r + d 
2 a 

r - d 
2 a 


1 - cos 


= 1 - cos 


Ei — E I 


+ arccos e cos 


Ei + E 


E 2 -E x 


+ arccos e cos 


2 

E 2 + E x 


(17) 

(18) 


In particular, for a real elliptical orbit to be possible the given scalar quantities 
must satisfy the relations r > d and 4a - r > d. 

If we define 


/?o = 2 arcsin 
then, using relation 


r + d 
4 a 


, y 0 = 2 arcsin 


r — d 
4 a 


1 - cos 9-2 sin' 


9 

2 ’ 


9 6 


and setting 


E 2 -E x ( E 2 + E i 

p =-r-+ arccos I e cos 


r = 


2 

Ei - E\ 


+ arccos e cos 


2 

Ei + E] 


(19) 

( 20 ) 

( 21 ) 


we find that the pairs 

(A 7) = 0ffo.ro). 03o,-yo), (2 n-/3 0 ,-y 0 ), (2n - (3 0 , y 0 ) 
satisfy equations (17), (18). 

Up to addition of the same integer multiple of 2n to both f5 and y, the pairs (21) 
are the only ones fulfilling (17), (18) and (16). 

From (19), (20) we obtain 

B + y Ei + E[ 

f3 - y = Ei - Ei, cos ——r— = ecos 


that yields 
In fact 


2 2 
nAt = (5 - y - (sinyS - sin y). 


n 0 • P-y P + y 

sin [3 - sin y = 2 sin —-— cos —— 


. Ei — E\ Ei + E\ . 

2e sin ——-cos ——-= e(sin E 2 - sin E x ). 
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□ 

The pairs (J3, y) given in (21) correspond to 4 geometrically distinct possible 
paths from the initial to the final position, see Figure 1. Given the points Pi, P 2 
and the attracting focus F, for a fixed value a of the semimajor axis, we find two 
different ellipses passing through P { and P 2 . They share the attracting focus F, 
but not the second focus (F* and F** in the figure). For each ellipse we have two 
possible arcs from Pi to P 2 , with different orientation, clockwise and counter¬ 
clockwise. 

The 4 cases are discussed in [6], [7], and are distinguished on the basis of the 
abscissa of the intercept Q of the straight line through P { and P 2 , on the axis 
passing through the foci of one of the ellipses, measured from its center. 



Figure 1: The four cases occurring in Lambert’s theorem. 

In [7] the 4 cases are also distinguished using the region K whose border is 
formed by the arc and the chord joining Pi and P 2 . We use this criterion for the 
classification given below. 

For a complete list of the equations coming from Lambert’s theorem, that need 
to be considered in our problem, we have to take into account the possible oc¬ 
curence of multiple revolutions along the orbit. Denoting by n the mean motion, 
the following expressions for At are obtained: 1 

i) At = Ti + Ikn/n, when the arc covers k revolutions and K contains neither 
of the foci; 

ii) At = T 2 +2kn/n, when the arc covers k revolutions, *R contains the attracting 
focus F but not the other one; 

'Here the region R is defined ignoring multiple revolutions. 
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iii) At = —T\ + 2 (k + 1 )n/n, when the arc covers k revolutions and K contains 
both foci; 

iv) At = —T 2 + 2 (k + 1 )n/n, when the arc covers k revolutions, K does not 
contain the attracting focus F but contains the other one; 

where T\,T 2 are given by 


nT 1 = A) - To - (sin Ad - sin y 0 ), 
nT 2 = A) + To - (sin A) + siny 0 ). 


The four cases above can be summarized in the equation 

77At = /3 - y - (sin /3 - sin y) + 2 kn, k e N, 
where the angles (3, y are defined by 


,i/3 r + d . 7 y r - d 
sin' - = ——, sin - 


2 4 a 


2 4 a ’ 


and 


0 < (3 - y < 2 7T. 

We also observe that in [8] there is a geometrical interpretation for the angles (3, 
7- 


12.2 Corrections to the observations 

We describe a procedure that could be used to correct the angular positions of a 
track by a pure geometrical argument. 

Assume we have the geocentric position vectors r j = p^f. + q j, j = 1... 4, with q 7 
the geocentric positions of the observer, from the radar observations of the celes¬ 
tial body. The vectors r 7 would be coplanar, if the orbit were perfectly Keplerian. 
In general this holds only approximately, due to the observational errors and to 
the perturbations which should be added to Kepler’s motion. We wish to correct 
these position vectors and define coplanar vectors r'., which are slightly different 
from Vj and keep the measured value pj of the topocentric radial distances. 

In the attempt to define a good approximation of the plane of this idealized 
Kepler motion, we compute the minimum of the function 

4 

v i-» Q(v) = Y/Xj ■ v) 2 

7=1 
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( 22 ) 


with the constraint |v| = 1. We obtain the equation 

4 

• yyrj - = », 

7=1 

with the Lagrange multiplier A e R, and consider the solution v min of (22) relative 
to the minimum eigenvalue A min . We take e v = v min as the direction of the Kepler 
motion plane, denoted by fl y . 

Then, for each j = 1... 4, we rotate the vectors pj = pp?. into a vector Kpj as 
follows (see Figure 2). 



Figure 2: Sketch of the correction of the line of sight. Here we skip the index j. 

Since we want to minimize the change in the line of sight, i.e. the observation 
direction e 1 ’, we rotate the latter around the axis orthogonal to the plane generated 
by e v , ej to reach the plane n v . In this way we draw a geodesic arc on the sphere 
with radius pj, centered at the observer position defined by q ; . This arc joins the 
position of the observed body with the plane n v . 

To describe this procedure in coordinates we introduce the angles Qj, <pj e 
[0,7r] defined by 

e v q; . v o 

cos6b =-, cos0; = e -e.. 

4j J 

The rotated vector %pj can be expressed as the linear combination 

%Pj = Aje v + B/ p 
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with Bj > 0 (since we do not want to rotate the line of sight by more than 90 
degrees). 

Now set the following conditions: 

(i) \Rpj\ = pj, 

(ii) [q 7 + K Pi ] ■ e v = 0. 

We obtain 

Aj + Bj + 2AjBj cos (pj = pj, 
qj cos 0 j + Aj + Bj cos (pj = 0. 

From the second equation we obtain 

Aj = -qj cos 6 j - Bj cos (pj, 

that substituted into the first yields 

Bj = — [ — Jpj - qj cos 2 6 j, 
sin (pj y ■> J 

so that 

Aj = - jj cos Qj + cot (pj ^pj - qj cos 2 d 7 j . 

We observe that this procedure works provided 

Pj > qj cos Qj, for j = 1... 4. 
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